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Abstract 

Thermally  conductive  composite  adhesives  are  desirable  in  many  industrial  ap¬ 
plications,  including  computers,  microelectronics,  machinery  and  appliances.  These 
composite  adhesives  are  formed  when  a  filler  particle  of  high  conductivity  is  added  to 
a  base  adhesive.  Typically,  adhesives  are  poor  thermal  conductors.  A  thorough  under¬ 
standing  of  heat  transfer  through  a  composite  adhesive  would  aid  in  the  design  of  an 
efficient  thermally  conductive  composite  adhesive. 

In  this  work,  we  provide  theoretical  foundations  for  use  in  design  of  thermally 
conductive  composite  adhesives.  For  proof  of  concept,  we  consider  a  two  dimensional 
model. 

We  prove  existence,  uniqueness  and  continuous  dependence  theorems  for  the  model. 

We  formulate  a  probability  based  parameter  estimation  problem  and  present  numerical 
results. 

Motivated  by  the  results  of  the  parameter  estimation  problem,  we  are  led  to  derive 
sensitivity  equations  for  our  system.  We  investigate  the  sensitivity  of  composite  sili¬ 
cones  with  respect  to  the  thermal  conductivity  of  both  the  base  silicone  polymer  and 
the  filler  particles.  Numerical  results  of  this  investigation  are  also  presented. 

Keywords:  thermal  conductivity,  composite  adhesives,  well-posedness,  inverse  prob¬ 
lems,  sensitivity  equations 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2001 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2001  to  00-00-2001 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Analysis  of  Thermal  Conductivity  in  Composite  Adhesives 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation, Raleigh, NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

49 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


1  Introduction  and  Motivation 


Adhesives  such  as  epoxies,  gels,  and  greases  have  numerous  commercial  and  industrial  ap¬ 
plications.  They  are  found  in  computers,  machinery,  home  appliances,  etc.  In  general,  these 
adhesives  are  very  poor  conductors  while  in  many  applications  it  would  be  advantageous  for 
them  to  possess  significant  thermal  conductivity.  Consequently,  researchers  have  been  study¬ 
ing  thermally  conductive  composites  or  filled  materials:  base  materials  such  as  epoxies,  gels, 
and  greases,  which  are  filled  with  thermally  conductive  particles.  Filler  particles,  such  as 
diamond  dust,  carbon  fibers,  or  aluminum  particles,  with  higher  thermal  conductivities  are 
added  to  create  a  composite  material  that  is  a  better  thermal  conductor  than  the  original 
material.  These  thermally  conductive  composites  could  then  replace  the  poorly  conduct¬ 
ing  adhesives  currently  in  use  in  applications  such  as  microelectronics,  circuit  boards,  heat 
exchangers,  machinery,  and  appliances. 

Adding  particles  with  a  high  thermal  conductivity  has  not  had  as  significant  an  impact 
on  the  overall  or  effective  thermal  conductivity  of  the  composite  as  anticipated.  In  order  to 
address  this  issue,  we  investigate  design  methodologies  for  these  composite  materials.  Our 
goal  is  an  improved  understanding  of  heat  transfer  through  a  composite  material  and  an 
increased  knowledge  of  the  impact  the  composite  design  has  on  this  thermal  process. 

The  goal  in  creating  a  thermally  conductive  composite  is  a  significant  increase  in  the 
thermal  conductivity  of  the  composite  over  the  thermal  conductivity  of  the  unfilled  material. 
There  are  several  design  considerations,  including  the  choice  of  particle,  the  particle  geometry 
and  the  size  and  shape  of  the  particles.  We  concentrate  our  presentation  here  on  a  composite 
material  with  a  fixed  geometry  and  consider  the  role  of  the  particles.  (For  the  effects  of 
varying  the  geometry  see  [8].) 

In  the  sections  below  we  present  a  mathematical  model  to  describe  the  heat  transfer 
through  a  composite  silicone.  We  show  that  the  mathematical  model  is  well-posed.  In  par¬ 
ticular,  we  show  there  exist  unique  weak  solutions  to  this  mathematical  model.  Furthermore, 
these  solutions  are  continuously  dependent  on  the  initial  conditions,  forcing  function,  and 
parameters. 


1 


We  formulate  a  probability  based  parameter  estimation  problem  based  on  results  in  [1] 
where  the  parameters  are  viewed  as  realizations  of  random  variables.  This  approach  allows 
for  uncertainty  in  the  model  parameters  as  well  as  the  data,  e.g.,  see  [9,  12].  We  introduce 
a  formulation  for  this  approach  in  the  context  of  our  model  and  present  some  numerical 
results. 

Finally,  we  rigorously  derive  sensitivity  equations  for  our  mathematical  model.  We  then 
numerically  solve  these  equations  and  show  that  the  model  is  more  sensitive  to  some  model 
parameters  than  others.  These  results  provide  insight  into  the  results  of  our  parameter 
estimation  problem. 

The  silicone  system  used  as  the  base  for  our  composite  silicone  has  a  thermal  conductivity 
of  approximately  0.12  W/mK  (Watts  per  meter-Kelvin).  The  base  silicone  consists  of  a  vinyl- 
functional  siloxane  (commonly  referred  to  as  a  resin  or  polymer)  and  a  hydride-functional 
siloxane  (commonly  called  a  crosslinker).  When  these  two  liquid  components  are  cured  the 
hydride  adds  to  the  double  bond  in  the  vinyl  group  to  form  a  linkage 

SiCH  =  CH2  +  H-Si - >  SiCH2-CH2Si 

which  is  sufficient  to  form  a  solid.  For  ease  of  reference,  we  refer  to  the  silicone  system 
as  the  silicone  polymer.  In  hopes  of  creating  a  composite  silicone  with  a  higher  thermal 
conductivity,  filler  particles  with  a  greater  thermal  conductivity  are  added  to  the  silicone 
polymer.  A  wide  variety  of  filler  particles,  including  aluminum  particles,  carbon  fibers  and 
diamond  dust,  can  be  added  in  varying  concentrations.  For  our  sample  composite  silicone, 
we  use  Grade  6  aluminum  which  has  a  thermal  conductivity  of  217  W/mK  and  concentrate 
here  on  composites  with  25%  by  volume  concentration  of  particles. 

There  are  a  variety  of  methods  available  to  measure  thermal  conductivity.  For  our  data 
collection  we  employed  a  Holometrix  Model  Microflash.  The  Microflash  uses  a  laser  flash 
method  which  allows  measurements  to  be  taken  at  room  temperature.  The  software  used  in 
conjunction  with  the  machine  is  Microflash-RT,  version  2.25. 

This  method  works  well  for  materials  of  uniform  density,  i.e. ,  non-composites.  However, 
we  are  using  composite  materials  and  it  is  known  that  this  results  only  in  some  measure  of 
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the  “average”  or  “effective”  thermal  properties  of  the  composite  material.  It  is  difficult  to 
answer  design  questions  about  the  composite  silicone  based  only  on  these  effective  properties 
of  the  composite. 

Analysis  on  a  single  test  piece  yields  the  “effective”  diffusivity,  specific  heat  and  thermal 
conductivity  of  the  sample  based  on  an  average  of  three  trials.  In  addition  to  these  three 
averaged  values,  the  Microflash  outputs  the  diffusivity,  specific  heat  and  thermal  conductivity 
of  each  trial  and  the  voltage  at  eight  different  times  for  each  of  the  three  trials.  The  eight 
times  recorded  are  the  time  in  milliseconds  (msecs)  to  reach  0,  20,  30,  40,  50,  70,  80,  and 
100  percent  of  the  temperature  rise  (equivalently  voltage  rise,  see  [8]).  For  further  details  on 
the  data  collection  method  and  experimental  results,  see  [8]. 

2  Problem  Formulation 

2.1  Model 

Since  the  fundamental  process  of  our  problem  is  heat  transfer  through  the  composite  silicone, 
the  foundation  of  our  model  is  the  transient  heat  equation  [11],  While  keeping  the  compo¬ 
sition  of  the  composite  silicone  and  the  data  collection  process  in  mind,  it  is  necessary  to 
make  a  few  simplifying  assumptions.  First,  we  assume  all  heat  from  the  heat  source  (a  laser) 
flows  through  the  composite  silicone  and  into  the  heat  sink  (an  IR  detector),  as  depicted  in 
Figure  1.  Second,  since  the  composite  silicone  slice  is  very  thin,  we  assume  there  is  no  heat 
loss  through  the  sides.  Thus  in  our  model  we  assume  the  sides  of  the  composite  silicone  are 
insulated.  We  use  a  flux  boundary  condition  to  describe  the  heating  on  the  source  side  of 
the  composite  silicone  due  to  the  laser.  On  the  sink  side  of  the  composite  silicone,  where 
the  IR  detector  is  located,  we  use  Newton  cooling  to  describe  the  boundary  condition  since 
that  face  of  the  composite  silicone  is  in  contact  with  the  ambient  air. 

However,  for  our  initial  study  of  the  problem  we  elected  to  reduce  the  three  dimensional 
model  to  a  two  dimensional  model  (solely  to  facilitate  numerous  computational  simulations). 
The  two  dimensional  model  can  be  thought  of  as  a  very  thin  interior  slice  (in  the  direction 
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Figure  1:  Three  dimensional  heat  transfer  model 


of  the  heat  flow)  of  the  three  dimensional  model  as  depicted  in  Figure  2.  We  assume  the 
composite  silicone  is  significantly  thicker  in  the  direction  normal  to  the  slice  compared  to  the 
slice  itself,  so  all  heat  will  flow  directly  through  the  composite  silicone  with  negligible  lateral 
dissipation.  For  our  experimental  test  pieces,  the  diameters  of  the  pieces  were  much  greater 
than  the  thickness  of  the  piece,  so  our  assumption  is  reasonable.  Furthermore  we  assume 
the  geometry  of  the  composite  silicone  is  uniform  normal  to  the  slice,  provided  we  are  in  the 
center  of  the  composite  silicone  away  from  lateral  edges.  Thus  the  two  dimensional  model 
heat  transport  properties  should  closely  resemble  that  of  the  three  dimensional  model. 
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Figure  2:  Two  dimensional  heat  transfer  model 

We  will  continue  to  assume  the  sides  of  the  composite  silicone  are  insulated.  The  bound¬ 
ary  at  the  heat  source  will  be  a  Neumann  boundary  condition  given  by  the  heat  flux  due  to 
the  laser  and  the  boundary  at  the  heat  sink  will  still  be  described  by  Newton  cooling. 
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It  is  clear  that  the  filler  particles  will  not  be  uniform  in  size,  and  will  most  likely  be 
randomly  dispersed  throughout  the  base  silicone  polymer.  However,  it  will  be  necessary  to 
know  the  size  of  each  particle  and  the  arrangement  of  the  particles  in  order  to  determine 
the  value  of  k ,  p,  and  cp  at  a  particular  point  in  the  composite  silicone.  To  facilitate  our 
modeling,  we  will  assume  the  filler  particles  are  fixed  and  comprise  the  appropriate  volume 
percent  of  the  composite  silicone  and  that  there  is  a  known  particle  arrangement. 

We  will  denote  the  ambient  temperature  by  7^  and  the  initial  temperature  of  the  com¬ 
posite  silicone  by  u0.  We  denote  the  Newton  cooling  constant  by  h.  and  define  S0(t)  to  be 
the  flux  due  to  the  heat  source.  Thus  if  u(t,  z)  is  the  composite  silicone  temperature  at  a 
given  time  t  and  coordinate  2,  we  have  the  following  system  describing  the  temperature  in 
the  sample: 

p(z)cp(z)u(t,  z)  =  V  •  (k(z)Vu(t,  z)),  z  e  Q 
k{z)j£{t,z)  |74  =  S0{t)  ( source ) 

<  ^)£(M|72  =  -  u(trz))\l2  {sink) 

k{z)^{t,z)  |7i  =  0 
k{z)^{t,z)  |73  =  0 
u( 0,  z)  =  Y( z),  2  6  fl 

s. 

where  O  =  [—  y,  y]  x  [— y,  y]  and  t.  G  [0,  T],  and  u  =  -|y,  with  ci,  c2,  and  T  assumed  finite, 
positive  constants.  Let  dkl  =  71  U  72  U  73  U  74  where 

-)»(»)  =  {(«,f):«€[-|,|]}. 

7 aM  =  {(- j,»)  :  s  e  [-|,|]},  and 

as  depicted  in  Figure  3.  We  define  Qs  C  Q  to  be  the  region  occupied  by  the  silicone  and 
C  O  to  be  the  region  occupied  by  the  filler  particles.  Note  O,  n  Op  =  0  and  Os  U  Op  =  O. 
It  is  important  to  note  that  (>.  cp  and  k  are  all  spatially  dependent.  They  will  have  one 
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Figure  3:  Two  dimensional  composite  silicone  slice 

value  in  the  silicone  polymer  and  another  value  in  the  filler  particles.  For  the  two  dimensional 
case  investigated  computationally  below  and  in  [8],  we  assumed  all  particles  were  circles  in 
a  known,  although  not  necessarily  uniform,  particle  arrangement. 

2.2  Matlab  PDE  Toolbox  Solutions 

The  major  benefit  of  using  the  two  dimensional  model  is  that  we  can  use  Matlab’s  Partial 
Differential  Equation  Toolbox  (PDE  Toolbox)  to  solve  (1).  Matlab’s  PDE  Toolbox  can  solve 
two  dimensional  parabolic  partial  differential  equations.  In  order  to  solve  (1)  using  the  PDE 
Toolbox  we  must  provide  the  boundary  conditions,  the  PDE  coefficients  and  the  composite 
silicone  geometry.  The  PDE  Toolbox  generates  a  triangular  mesh  using  the  Delaunay  tri- 
angulation  algorithm  and  numerically  solves  the  PDE  using  the  finite  element  method  with 
linear  elements  (the  only  type  of  elements  the  PDE  Toolbox  employs).  The  PDE  Toolbox 
automatically  defines  the  mesh,  although  the  user  has  the  option  to  refine  the  mesh.  See 
[14]  for  further  information  about  Matlab’s  Partial  Differential  Equation  Toolbox. 

Matlab’s  PDE  Toolbox  allowed  us  to  carry  out  simulations  for  many  different  geometry 
configurations.  For  example,  geometries  with  uniform,  shifted  and  random  geometries  can 
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be  considered.  We  can  (and  did)  also  consider  geometries  of  both  same  sized  and  varying 
sized  particles.  A  summary  of  our  numerical  simulations  for  different  geometries  (random, 
fixed,  etc.)  and  different  distributions  of  particle  size  is  given  in  [8]. 

For  the  numerical  results  we  present  in  this  paper,  we  assume  all  aluminum  filler  particles 
in  the  two  dimensional  model  are  circles  of  uniform  diameter  arranged  in  a  uniform  geometry, 
i.e. ,  the  particles  are  uniformly  spaced  and  aligned  in  rows.  Since  our  sample  composite 
silicone  contains  Grade  6  aluminum  particles  we  will  use  the  mean  diameter  of  the  volume 
distribution  provided  by  the  aluminum  supplier,  24.14  //  (microns),  as  the  diameter  of  each 
particle.  We  will  concentrate  on  the  25%  by  volume  composite  silicone,  but  all  ideas  presented 
here  extend  in  a  natural  way  to  composite  silicones  with  different  compositions. 

The  silicone  polymer  used  as  the  base  for  the  composite  silicone  wets  well,  meaning  it 
forms  a  thin  film  around  each  of  the  particles.  The  film  formed  around  each  particle  is 
estimated  to  be  approximately  50  angstroms.  Hence,  we  assume  each  particle  is  separated 
by  a  minimum  distance  of  0.01  microns  and  that  no  particles  touch  the  boundary  of  the 
composite  silicone  slice  which  is  321.5  //  wide  and  1638  //  high,  i.e.,  c\  =  321.5  and  c-i  =  1638 
in  the  two  dimensional  model  of  the  previous  section.  The  composite  silicone  slice  we  use 
in  computations  reported  on  below  contains  288  circles  of  diameter  24.14  ^ i  in  a  uniform 
arrangement  representing  the  25%  by  area  aluminum  particles  used  by  the  PDE  Toolbox  as 
the  geometry  for  our  composite  silicone. 

We  note  that  k,  p  and  cp  are  all  spatially  dependent.  In  order  to  differentiate  between 
the  value  of  each  parameter  in  the  silicone  versus  the  aluminum  particle  we  will  quantify 
this  variability  as  follows:  for  z  £  O,  the  value  for  each  parameter  is  given  by: 

z  £  Os 

Z  £  Op 

z  £  Os 

Z  £  Op 


k{z)  = 


p{z)  = 


Ps 

Pp 
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and 


I  cPs  z  E  fis 
cp(z)  =  ^ 

cpp  z  £  Qp  ■) 

where  ks,  kp ,  ps,  pp ,  cPs ,  and  cPp  are  all  finite  constants.  Observe  that  k,  ft.  and  cp  are  each 
functions  from  O  to  M  and  each  is  piecewise  constant. 

We  assume  the  composite  silicone  was  initially  at  the  ambient  temperature  and  the 
temperature  was  uniform  throughout  the  sample.  Thus  we  set  T(^)  =  7%,.  We  choose 
the  Newton  constant  to  represent  air  cooling  as  we  have  in  our  model.  The  exact  model 
parameters  used  in  the  simulations  reported  here  are  in  Table  1.  In  the  table,  g/cm3  is 
grams  per  cubic  centimeter  and  J/gK  is  Joules  per  gram-Kelvin. 


ka 

0.12  W/mK 

kp 

217  W/mK 

Ps 

1  g/cm3 

Pp 

2.7  g/cm3 

CPs 

1.55  J/gK 

CPp 

0.90  J/gK 

T 

oo 

296.15  K 

h 

350 

St 

4.32  x  107  W/m2 

Table  1:  Model  parameters 

The  source  flux  will  approximate  the  energy  in  the  laser  pulse.  The  laser  energy  as 
configured  in  the  Microflash  is  approximately  7  .J.  For  our  testing  there  is  a  20%  filter 
screen,  so  the  actual  laser  energy  is  approximately  1.4  J.  In  addition,  since  the  pieces  are 
graphite  coated,  technicians  at  Holometrix  estimate  there  is  an  additional  20%  energy  loss. 
Given  that  the  length  of  the  laser  pulse  is  330  microseconds  and  the  diameter  of  the  laser  is 
10  mm,  we  calculate  the  flux  due  to  the  laser  pulse  to  be  S(  =  4.32  x  10'  W/m2  (Watts  per 
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meter  squared).  Thus  the  source  flux  is  given  by 


S0(t) 


Se  0  <  t  <  tp 
0  tp  <  t. 


where  tp  =  0.000330  seconds.  While  we  believe  this  is  a  high  estimate,  it  is  sufficient  for  our 
purposes  here  and  in  [8]. 


3  Well-Posedness 

3.1  Problem  in  Variational  Form 

In  this  section  we  investigate  theoretical  issues  relating  to  the  two  dimensional  model  for¬ 
mulated  in  Section  2.  We  define  a  class  of  abstract  parabolic  equations  and  establish  that 
this  class  of  equations  is  well-posed.  Furthermore,  we  verify  that  our  two  dimensional  heat 
transfer  model  fits  this  class  of  equations  for  a  large  number  of  examples  with  different 
particle  shapes,  size  distributions  and  geometry.  These  results  guarantee  the  existence  and 
uniqueness  of  weak  solutions  to  our  model  as  well  as  continuous  dependence  on  the  initial 
data,  forcing  function  and  model  parameters  for  most  examples  of  interest.  We  note  here 
that  these  theoretical  results  (and  those  of  Sections  4  and  5)  are  easily  extended  to  three 
dimensions ,  but  we  choose  to  concentrate  on  the  two  dimensional  problem  since  it  relates  to 
our  use  of  Matlab’s  PDE  Toolbox. 

3.1.1  Preliminaries 

We  begin  with  the  two  dimensional  model  (1)  from  Section  2  with  rather  general  but  fixed 
particle  shapes,  sizes  and  geometry  of  location.  We  define  for  simplicity  of  notation  g(z)  = 
p(z)cp(z)  throughout  and  assume  there  exists  constants  Rl  and  Rjj  such  that 

0<RL<  g{z)  <  Rjj  <  oo  (2) 

and  constants  KL  and  KtJ  such  that 

0  <  Kl  <  k(z)  <  Kj  <  oo  (3) 
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for  all  z  G  Q.  Furthermore  we  assume  h  is  a  finite,  positive  constant. 

Define  H  =  L2(Q)  with  the  usual  L 2  inner  product,  (•,  ■)i2,  and  define  IK  =  L2(Q)  with 
the  weighted  inner  product  =  (g  Note  the  norms  generated  by  the  fiC-inner 

product  and  the  L2-inner  product  are  equivalent.  Let  V  =  iL1(D)  with  inner  product 

(</»,  i’)v  =  <V0,  V'0)L2  +  (0,  0)l2 

for  0,0  G  D  and  let  V  =  iL1(0)  with  inner  product 

(0, 0)v  =  (V0,  Vr):l|  +  (0, 0)L2  (4) 


for  0,  -0  G  V.  We  will  want  to  use  the  following  equivalent  representation  of  the  inner  product 
in  V: 


(0,0)?  = 


(V0,  V0K  + 


(Tr2  (j)){z){Tr2  i>)(z )  dS2 


(5) 


where  Tr,-  :  iL1(D)  -G  (yf)  is  the  continuous  trace  operator  mapping  /  G  HL(Ll)  -G 
C  iL°(7j)  =  on  7,-,  j  =  1,2,  3,  4,  with  (Trj  f){z)  =  f{z) |7..  The  following 

theorem  from  Maz’ja  [16,  p.  27]  shows  the  norms  generated  by  these  inner  products  are 
equivalent: 


Theorem  3.1  (Maz’ja)  Let  Q  be  a  bounded  domain  in  Rn  such  that  L{p(Q)  C  Lp(Q).  Let 
30'u)  be  a  continuous  functional  in  Wp(Q),  ^(n^i)  ^  0  for  any  nonzero  polynomial  Iii-\  of 
degree  not  higher  than  l—  1.  Then  the  norm 


||Vyu||£p(Q)  +  T(u) 

is  equivalent  to  the  norm  in  Wp(Q). 

Here  Lp(il)  is  the  space  of  distributions  on  Q  with  derivatives  of  order  t  in  LP(Q),  and 
Wp(Q)  =  Lp(Q)  fl  Lp(Q).  Also  =  { Da },  where  |o|  =  l  for  a  a  multi-index  (au,  •  •  •  ,an) 
with  Da  =  Djj1  •  •  •  Df™ .  (Note  for  our  problem  11': ( 0 )  =  110 (Q)  =  iL1(D),  i.e.,  p  =  2  and 
1=1  and  n  =  2.)  If  we  define 

T('u)  =  y  |£r(u)|2cLc 
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where  I  c  911  and  tr  is  the  continuous  trace  operator  mapping  u  G  H 1(H)  — »  H 2  (911)  C 
If0 (9 11)  =  L2(911),  then  we  have 

|9r(«)|  =  y  |tr (w) |2r/rc 

<  /  |lr(u)|2c/.x 

J  an 

<  Ii  ||«||#1(Q) 

and  hence  T  is  a  continuous  functional  on  H"1(H).  Also  note  that  taking  t  =  1,  J7_4  =  c, 
where  c  is  any  non-zero  constant,  and 

T(e)  =  y  | tr  (c)  1 2cZrr 

is  non-zero  if  and  only  if  the  measure  of  T  is  non-zero.  Thus  the  norms  generated  by  the 
inner  products  (4),  (5)  are  equivalent  by  Theorem  3.1  and  there  are  finite,  positive  constants 
Ml,  Mu  such  that 

mlm\1<  m\%  <MuM\l  (6) 

for  all  0  G  V  where 

M\l=\\Vct>\\ [  |Tr2  0|2  clS2,  (7) 

j  72 

is  the  norm  generated  by  (5). 

Define  a  sesquilinear  form  a  :  V  x  V  — >■  M  by 

a(0,  0)  =  <-£V0,  V0)jc  +  h  [  (Tr2  0)(7)(Tr2  0)(7  ^2.  (8) 

9  J  72 

Note  a  defines  an  operator  A  G  Jzf(V,  V*)  where  (A0,  0)v*,v  =  <7(0,0)  and  Jz?(V,  V*)  is  the 
set  of  all  bounded,  linear  functionals  from  V  to  V*.  This  follows  due  to  the  continuity  of  a 
on  V  x  V  guaranteed  by  (6)  (see  (12)  below). 

Define  F  :  [0,T]  — >■  V*  by 

[F(1)](0)  =  hT^  I  (Tr2  0)(,z)  dS2  +  50(t)  /  (Tr4  0)(F)  dS4  (9) 

7  72  I74 

for  (/>  G  V. 
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3.1.2  Weak  Solution 


Suppose  u  solves 


i.e.,  for  all  D  £  V, 


By  definition  we  have 


u  +  Au  =  F  in  V* , 

{u  +  Au  -  F.  l[’)V*y  =  0. 


(10) 


{u{t),  V’)v*,v  —  (— Au(t ),  V’)v*,v  +  (U  V')v*,v 


=  < — kVu(t),  W>)rK  -  h  /  (Tr2  w(t))(z)(Tr2  ^)(  :)  rIS'2 

5  .7  72 

+  hTco  f  (Tr2  *>)(:)  r/.S2  |  .%(/)  f  (Tr4  4>)(z)  dS4 

j  72  774 


=  {—kVu(t),  V0)L2  -  /,  /  ((Tr2  «(*))(*)  -  T00)(Tr2  </’)(*)  dS2 

<7  72 

+  So(t)  /  (Tr4  V’)(-)  dS4- 

J74 

Now,  if  u  G  y  and  A:  V«  G  V,  using  the  Divergence  Theorem  and  the  vector  identity 


V  •  (k(z)Vs(t,z)4>(z))  =  (V  •  (k(z)Vs(t,z))4>(z)  +  k(z)(Vs(t,  z)  ■  Vt/’(s)) 

one  can  argue  that  a  solution  u  of  (10)  (if  it  exists)  is  a  weak  solution  of  (1).  That  is,  (10) 
is  the  weak  form  of  (1). 


3.2  Well-Posedness  (Existence,  Uniqueness,  Continuous  Depen¬ 
dence  on  Data) 

We  establish  existence  of  solutions  to  parabolic  systems  of  the  form 


u  +  Au  =  F  in  V* 
u(0)  =  u0. 


(11) 
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We  will  use  the  Gelfand  triple  V  'K  =  fit*  V*  where  the  embedding  V  IK  is  dense 
and  continuous  with  ||0||j{  <  c||0||v  for  all  0  G  V  and  some  finite  constant  c  >  0.  Note  the 
desired  inequality  holds: 

|.[0||v  =  (V0,  V</»)rK  +  (0,  4>)l* 

>  (0,0)  L2 

=  (  0i  0)jC 

9 

—  R  (  /  1 1 0 1 1  ac 

so  in  fact  ||0||w  <  VSfo  1 10|  |v- 

We  shall  argue  and  then  use  two  standard  conditions  on  a  defined  by  (8): 

(1)  The  form  a  is  V-bounded:  for  all  0,0  G  V,  there  exists  a  B  <  oo  such  that 

|a(0,0)|<B||0||v||0||v.  (12) 

(2)  The  form  a  is  V-coercive:  for  all  0  G  V,  there  exists  a  C  >  0  such  that 

M0,0)|  >  C'||0]|v-  (13) 

We  also  find  that  if  the  source  So  is  in  L2(0,T),  then  the  forcing  term  F  defined  by  (9) 
satisfies 

F  G  L2(0,  X;  V*).  (14) 

To  show  that  a  of  (8)  satisfies  (12),  we  use  (2),  (3),  and  (7): 

M0,  0)|  =  |<-fcV0,  V0)^  +  /;  /  (Tr2  0)(z)(IV2  4>)(z)  dS2\ 

9  F/2 

<  max{^1/W,/i}(||V0||J{||V0|k+  ||Tr2  0||72||Tr2  0||72) 

<  max{^1AV,/i}(||V0||^+  ||Tr2  0||7J(||V0||^  +  ||Tr2  0||72) 

<  max{JR“1AV,/i}(2||0||?)(2||0]|?) 
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<  4 Ml  ma x{Rl  Ki;.  h) \ \<p\  |v|  |V’|  |v 


where  \\f\\l2  =  /  \f{z)\2  dS2  for  any  /  G  H 2(72).  Thus  o  is  V-bounded. 


'72 


Similarly,  for  (13)  we  have  the  following: 


a(0,  </>)  =  /-IVo.  V<j>)oi  +  /i  /  (Tr2  0)(^)(Tr2  </>)(^)  dS2 
9  J  72 


>min{JR-iA'i,/r}(||V0||^+  /  |(Tr2  </>)(T)|2  dS2) 

J  72 


>  /'MV  0  lv 


so  a  is  V-coercive. 

In  order  to  see  that  (14)  holds,  recall  Tr2  :  V  =  H 1(fi)  — >■  H^{^2).  So,  for  any  (/’  G  V, 
Tr2  F  E  H%( j2)  C  L2(j2),  and  hence  I  Tr2  </’  clS2  E  M.  Thus  i[>  — >■  J  Tr2  A  is 

•272  72 

a  continuous  mapping  from  V  — >■  M,  i.e.,  it  is  in  V*.  If  S0  E  L2(0,T)  (which  we  assume 
throughout),  then  F  E  L2(0,X;  V*). 

Given  the  above  hypothesis,  the  system  (11)  is  equivalently  written 


(«(*),  #  +  ''/’)  =  V’) 

u(0)  =  Mo 


(15) 


for  </>  G  V  where  the  duality  product  (•,  •)  is  (•,  -)v*,v- 

Assume  for  the  moment  that  (15)  has  a  solution  u.  We  derive  an  a  priori  bound.  Let 
'W  =  u(t )  for  a  fixed  t.  Substituting  into  (15)  we  obtain 


(u(t),u(t))V*y  +  <T{u{t),u{t))  =  (F(t),u(t))V*y 
and  since  (u(t),  u(t))v,,v  =  ||{||V)llV  we  see 

+  =  (F(f),u(i))v*,v 

for  any  t  in  a  given  interval  [0,  T].  Integrating  from  0  to  t  we  have 

£{~{"«(o"t} +*(«(«, «({))}  <-;« = <if 
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and  thus 


\\\um\l±  <*e  =  |V(8,«(f))w 

Using  (13),  the  Cauchy  Schwartz  inequality  and  the  fact  2 cib  <  a2  +  b2,  we  have 

IM*) Woe  +  2C  r  II «(Ollv  #  <  Ifwollw  +  2|  /V(0, «(0>v,v  del 

Jo  Jo 

<  H«o||^4  2  /V(Ollv*||u(Ollvd£ 

Jo 

<  iw&  +  ^Vk)HI.  de  +  c^'n«K)iiide 

and  hence 

umul  +  c  J‘ \\U{o\& as  <  \\'uo\\l<  +  ^  j‘ mewl.  <%.  (ie> 

Thus 

Mmli  +  c  f  IMOWI  dt  <  c 

Jo 

where  C  =  (7( | |w0|  | jc,  C,  |  |iH  Ii>(o,t;v*))- 

The  a  priori  bound  arguments  are  the  basis  of  existence  as  well  as  continuous  dependence. 
Using  them  along  with  quite  standard  arguments  (see  Chapter  III  of  [13],  §26  of  [23]),  we 
can  establish  the  desired  existence  and  uniqueness  (detailed  arguments  are  given  in  [8]). 

Theorem  3.2  Under  assumptions  (12),  (13),  and  (14),  for  u0  G  V,  there  exists  a  solu¬ 
tion  of  (11)  (and  hence  a  weak  solution  of  (1))  with  u  G  L2(0,T;  V)  and  u  G  L2(0,T;  V*). 
Furthermore,  this  solution  is  unique. 


In  a  similar  standard  approach  (again  see  [8,  13,  23])  one  can  establish  continuous  de¬ 
pendence  of  solutions  on  initial  data  and  forcing  function.  We  only  recall  the  ideas  here. 
Suppose  n  //( • ;  uq,  F)  is  a  weak  solution  of 


u  +  Au  =  F 
U{  0)  =  U0 


(17) 


IB 


and  suppose  un 


un (  • ;  un o,  Fn)  is  a  weak  solution  of 


iin  -b  Aun  Fn 

Wn(0)  ^n0 


(18) 


with  uno  — >■  ?/o  in  IK  and  Fn  F  in  L2(0,T;  V*).  Given  systems  (17)  and  (18)  we  see 


(«  -  «n,  V’)  +  cr(w  -  V>)  =  (F{t)  -  Fn{t ),  V’) 

w(0)  -  un(0)  =  u0  -  un o 


(19) 


for  all  V’  €  V.  If  we  let  0  =  w(£)  —  un(t )  in  (19)  and  use  the  same  arguments  as  in  establishing 
(16),  we  see 

\\u{t)  -  Un{t)\\x  +  C  f  NO  -  «n(Ollv  <  \\U0  ~  Uno\\lt  +  ^  1 1^(0  “  |  |v* 

and  thus 

ft  i 

I|«(t)  —  Un{t)\\l{  +  C  J  ||(f)  —  un{  Ollv  G  1 1  «0  —  W'nO  1 1 (K  +  ^||-f  ~  Fn  |  |l2(0,T;V*)  • 

Thus  given  that  un o  — >■  «o  in  IK  and  Fn  — >■  F  in  L2(0,T;  V*),  we  see  u  — >■  un  in  (7(0,  T;  IK) 
and  also  in  L2(0,T;  V).  Indeed  we  have 


Theorem  3.3  The  mapping  (uo,F)  — >■  u(  • ;  uo,F),  where  u(  ■ ;  uo,F)  is  a  solution  to  (11) 
is  continuous  from  IK  x  L2(0,T;V*)  to  (7(0,  T;  IK)  fl  T2(0,  T;  V). 


L2 


We  remark  that  one  can  actually  establish  the  somewhat  stronger  continuity  from  IK  x 
(0,  T;  V*)  to  It  =  L2(0,  T;  V)  fl  iP(0,  T;  V*),  see  [13,  23]  for  details. 


3.3  Continuous  Dependence  on  Parameters 

3.3.1  Continuous  Dependence  on  k  and  h 

Suppose  u(  ■ ;  k,  h)  is  a  weak  solution  of 

{u  +  Au  =  F 
w(0)  =  u0 


(20) 
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and  let  un(  ■ ;  kn,  hn)  be  a  weak  solution  of 


Un  “b  A n  un  Fn 

Un  ( 0 )  7/  0 


(21) 


where 


(An<i>,  4’)v*,v  =  crn{4 >,  4 ’)  =  (~knV<j),  V0)jc  +  hn  I  (: Tr2  <j>){z)(Tr2  4>){z)  clS2  (22) 
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and 


(Fn{t),  4’)v*y  —  KToo 


{Tr2  b)(  :)  ri52  +  50(t)  /  (Tr4  /;)(c)  d54. 


'72 


'74 


Let  {//„}  be  a  sequence  such  that  hn  — )■  h  and  let  {/cn}  be  a  sequence  such  that  kn  — )■  /,; 
uniformly  in  2.  i.e.,  kn  — )■  /,;  in  C(O).  Since  we  know  hn  — )■  h,  this  implies  Fn  — )■  F 
in  L2(0,T]V*).  From  the  previous  section  we  know  solutions  depend  continuously  on  the 
forcing  function  when  Fn  — >■  F  in  L2(0,T;V*).  Thus  without  loss  of  generality  we  can 
suppress  the  dependence  of  F  on  hn  and  take  F  for  Fn  in  (21)  with  un(-',knJin)  a  weak 
solution  of 


un  T  Anun  F 

7'  ri  (11)  77  o 


(23) 


where  Tn  is  defined  as  in  (22).  By  (3)  we  know  for  Ni  sufficiently  large  there  exist  constants 
Kl  >  0  and  Kjj  <  oo  such  that  kn ( z)  G  [KL.  K/;].  k(z )  G  [Kl,Kjj]  for  all  2  G  O  and  all 
n  >  JV4.  Since  hn  — )■  //.  for  i¥2  sufficiently  large  there  exist  constants  hL  >  0  and  hy  <  oo 
such  that  hn  G  [hr,.  hu\,  h  G  [hr.  -  hu]  for  all  n  >  N2.  Without  loss  of  generality,  hereafter  we 
assume  all  n  will  satisfy  n  >  N  =  max{Ari  .Ay}. 

Note  an  is  uniformly  (in  n)  V-coercive  satisfying  rr„(o.  o)  >  C'i||</)||^,  where 


C,  =  mm{R-lKL.hj,}Mk] 


and  Ci  is  independent  of  n. 

Subtracting  the  weak  form  of  (23)  from  (20)  we  obtain 

(u{t)  -  un{t),  4>)v*y  +  cr{u{t),  </>)  -  an{un{t),  4’)  =  0 


17 


for  all  'ip  G  V.  Thus 


(u(t)  -  iin{t),  i’)v,v  +  {-(kVu{t)  -  knVun{t)),Vi’)oc 

( h(Tr2  u{t))(z)  -  hn{Tr2  un{t)){z)){Tr2  4>)(z)  dS2  =  0 


72 

for  all  A  G  V. 

Adding  and  subtracting  terms  we  see 

(u(t)  -  Un(t),  4>)vy  +  -  *n)Vu(t),  Vt/’):K  +  {^kn{V U,{t)  -  Vun(£)),  0}.<K 

+  (h  -  hn)  I  ( Tr2  u(t))(z)(Tr2  4>){z)  dS2 

J  72 


+  K  /  (Tr2  (u(t)  -  un(t)))(z)(Tr2  4>)(z)  dS2  =  0 

J  72 

for  all  4’  £  V. 

Fix  t  G  [0,T]  and  let  'W  =  u(t )  —  in  the  previous  equation.  Thus 

“{f|«(^)  -  «,»(*)!&}  +  (^MVu(t)  -  Vu„(i)),  Vti(t)  -  Vitn(i))M 


K  /  (Tr2  («(£)  -  un(t)))(z)(Tr2  (u(t)  -  un(t)))(z )  dS2 


'72 


=  (-(&n  -  k)Vu{t),  Vu{t)  ~  Vun(t));H 
9 

+  {hn  -  h )  /  (Tr2  u(t))(z)(Tr2  (u(t)  -  un{t)))(z)  dS2 

J  72 

and  hence  using  the  definition  of  crn  and  integrating  from  0  to  t,  for  t  G  [0,  T]  we  obtain 


IK*)  -  «»»(*)!&  +  2  /  <7„(u(f)  -  “  «n(f)) 

JO 

=  2  [  (-(kn  -  k)Vu(£),  Vu(f )  -  Vun(f  ))*; 

Jo  9 

+  2  /  (hn-/i)  [  (Tr2  u(i))(z)(Tr2  (u(0  ~  un(4)))(z)  dS2  d£. 

JO  J  72 

Using  the  fact  <rn  is  V-coercive  on  the  left  side  of  the  equation,  and  the  Cauchy-Schwartz 
inequality,  the  relation  2 cib  <  a2  +  b2,  and  the  equivalent  form  of  the  V-norm  generated  by 
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(5)  on  the  right  side,  we  see 

\\u{t)-un{t)\\l<  +  2Cl  I  ||u(f)  -  «n(Ollv  d€ 

Jo 

<  2  fw-iK  -  A:)V«(e)|k||V«(e)  -  V«„(Olk  cli 

Jo  9 

+  2  [  | (hn-h)  I  (Tr2u(0)(z)(Tr2(u(0-Mm(z)dS2\d£ 

Jo  J  72 

<  fr  J‘ri2  tt*,  -  *ii?oIIv«(0||m  <k  +  y  J‘  nv«(f)  -  v«„(e)iiif  <k 

+  ~ttl~  /  l^-/f  /  |(IV2u(0)(^)|2d52de 

^'1  </0  J  72 

+  J‘  j  IPh  MO  -  «»(f)))WI*  df 

^  _ *11°°/  <;f  +  y/  H“K) _ «»(?)llv <(£ 

+  -  *i2 /'  worn  <k  +  yJ‘  H“(f)  -  “”(?)iiv  rff 

where  1 1  •  |  loo  is  the  C(Q)  norm.  Combining  like  terms  we  see 

\\u{t)  -  UnWWx  +  Ci  I  ||«(0  -  «n(f)||v 

Jo 

<  (jt|j life  -  k\\l  +  K  -  fc|»)  J‘  WOlll  <%■ 

If  we  define  Gn(T)  by 

Gn{T)  =  (ctfrWkn  -  kW2^  +  c  M2  I  hn  ~  ^|2)||«||l2(o,t;V) 

1  L  1  L 

and  recall  kn  —>  k  in  C(O),  hn  —>  h  and  u  G  L2(0,X;  V),  we  see  Gn(T)  — >■  0.  Thus 

I  W«)  -  «„WII«  +  y  f  MO  “  MOIII rf{  £  G”(T)  (24) 

and  so  as  n  — >■  oo,  u  — >■  un  in  C(0,  T;IK)  and  in  X2(0,T;V).  Thus  the  solution  depends 
continuously  on  the  parameters  A;  and  h.  We  actually  have  proved  the  somewhat  stronger 
result: 
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Theorem  3.4  The  mapping  (k,  h )  >  u(  ■  ■:  k ,  // ),  where  u(  ■ ;  ft,  /i)  is  a  weak  solution  to  (20) 
is  Lipschitz  continuous  from  C(Q)  x  to  G(0,T;!K)  D  L2(0,T;  V). 


3.3.2  Continuous  Dependence  on  p  and  cp 

As  before,  let  g  =  pcp  and  define  a  sequence  {gn}  such  that  gn  — >■  g  uniformly  in  z,  i.e. , 
in  C(Q).  Thus  we  know  by  (2)  for  N  sufficiently  large,  there  exist  constants  Rl  >  0  and 
Rjj  <  oo  such  that  for  N  sufficiently  large,  gn(z )  G  [R-r,-  Ru]-  g{z)  G  [RL,Ru]  holds  for  all 
n  >  N  and  all  zGO.  Without  loss  of  generality,  from  hereafter  we  assume  all  n  will  satisfy 
n  >  N. 

Let  TCn  =  iL1(f2)  with  the  weighted  inner  product  (•,  -)^n  =  (gn ■ ,  -)L2.  Note  the  J~Cn  norm 
is  equivalent  to  the  J~C  norm  uniformly  in  n  and  there  are  finite,  positive  constants  Ju 
such  that 

Wilt  <  \mn  <  JuM&  (25) 

for  all  n.  Define  Vn  =  H1( Q)  with  inner  product 

(0G0)v„  =  (V <x 


which  has  equivalent  representation 


(0,  f,)Vn  -  (V0,  VV>):k„  +  /  h(Tr2  c/>)(z)(Tr2  4>)(z)  dS2 

j  72 

by  Theorem  3.1. 

Suppose  u(- ;  g)  is  a  solution  of 

u  +  Au  =  F 

u{  0)  =  «0 

in  V*  and  suppose  un ( • ;  gn)  is  a  solution  of 

fin  T  An  un  Fn 
^n(O)  «0 


(26) 


(27) 


20 


in  V*  where 


(An(t),  -ij')v*yn  =  crn(<t>,4>)  =  (— kV(f>,  Vt/’)w„  +h  (Tr2  4>)(z)(Tr2  4>){z)  dS2 

9n  J  70 


{Fn(t)<j),  V’)v*,v„  —  ( F(t)(f),.ij))v*y . 

for  all  (/’  G  Vn  =  V. 

Subtracting  the  weak  form  of  (27)  from  (26)  we  obtain 

i’)v*,V-(un(t ),  #)v*,Vn  +  CT(w(t),  '(/’)  -  Vn{Un(t),  4’) 

=  (/•'(/).  ?l)v.v  (Fn{t),  4’)v*yn 

for  all  xp  G  V.  Thus 

(«(/). i>)y.v  -  («„(*),  V’)v*,v„  +  (^&Vu(i),  VV>)rK 

+  h  [  (Tr2  u(t))(z)(Tr2  4>)(z)  dS2  -  {-kWun(t), 

J 72  9n 

-  h  I  (Tr2  un(t))(z)(Tr2  4>){z)  dS2  =  0 

J  72 

and  hence 


(«(*), V’)v*,v  -  (hn(£),V’)v*,v„  +  (A:V«(t),  V'0)L 2 

+  /r  /  (Tr2  u(tp)(z)(Tr2  4)(z)  dS2  -  (kVun(t),V4’)L* 


-hi  (Tr2  un(t))(z)(Tr2  4>)(z)  dS2  =  0 


for  all  4’  G  V.  Define  <r  by 


a(</>,  0)  =  V r)/;.!  +  /;  /  (Tr2  </»)(^)(Tr2  </’)(*)  ^2 


and  note  a  is  V-eoercive  satisfying  <r(</>,  0)  >  <772 1 1 <^>|  | v  where 


C2  =  min{.Rj71i£i, 
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Adding  and  subtracting  terms  we  see 

(u(t),  V’)v*,v  -  (u(t),  i’)v*n,v^  +  (u(t)  ~  un(t),  '0)v*,v„  +  a{u(t)  -  un{t),4>)  =  0. 

Fix  t  G  [0 ,  T]  and  let  ip  =  «(£)  —  un(t)  in  the  previous  equation.  Thus 

7^{IM*)  -  Un{t)\\xn}  +  <*(«(*)  -  «n(*),  «(*)  -  «n(*)) 

=  (yU(i)  -  U{t),u{t)  ~  Un(t))V*y 

and  hence  integrating  from  0  to  t  we  obtain 

||'u(t)  -  u„(i)|&n  +  2  [  a(u(£)  -  un(£),u(£)  -  u„(f))  d£ 

Jo 

=  2  [  ((— -1  )u{£),u(t)-un{t))v*y  d£. 

Jo  9 

Since  the  J~Cn  norm  is  equivalent  to  the  J~C  norm  uniformly  in  n,  using  (25)  we  can  rewrite 
this  equation  as 

JL\\u(t)-un(t)\\li  +  2  I  a(u(£)  -  tt»(f),«(0  -  ttn(0)  d£ 

Jo 

<  2  f  |((—  -  l)u(^),u(t)  -  un{t))v*y\  cl£. 

Jo  9 

Using  the  fact  that  a  is  V-coercive  on  the  left  side  of  the  equation,  and  the  Cauchy  Schwartz 
inequality  and  the  relation  2 cib  <  a2  +  b 2  on  the  right  side,  we  see 

JL\\u{t)-un{t)\\lc  +  2C2  I  ||«(0  -  «n(Ollv  d £ 

Jo 

<2  r  in— -  i)«(f)Mviwo  -“»(f)Hv  rff 

Jo  9 

<  A  f  IK  V  -  D*(Ollv.  rff  +  C,  f  ||«(f)  -  Ii„(0lll  d(. 

c'2  .7o  9  Jo 

Combining  like  terms  we  see 

JL\\u{t)  -Un{t)\\l<  +  C2  f  ||u(f)  -  «n(C)Hv  dt,  <  TT  [  \\—  -  MlUHOWv*  d£. 

Jo  c'2  Jo  9 
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Let  Gn(T )  be  defined  by 


Gn(T) 


1  II  9n-  9 
G'2  9 


2 

oo 


w|U2(0,T;V*) 


and  note  Gn(T )  — >■  0  since  gn  g  in  C(il)  and  u  G  L2(0,T;  V*).  Hence  we  have 

JL\\u(t)  ~un(t) Wx +  C2  I  ||u(f)  -  «n(Ollv  <  Gn{T) 

Jo 

and  so  un  — >■  u  in  C(0,T;IK)  and  in  L2(0,T;V).  Thus  the  solution  depends  Lipschitz 
continuously  on  the  parameters  p  and  cp.  We  have  proved  the  following  theorem: 


Theorem  3.5  The  mapping  g  — >■  u(-;g),  where  u(-\g )  is  a  solution  to  (26)  is  Lipschitz 
continuous  from  C (fi)  to  C(0,  T;  TC)  fl  L2(0,T;  V). 


4  Formulation  of  the  Inverse  Problem 

4.1  Preliminaries 

For  {/,  }''  .  C  [0,T],  T  <  oo,  we  can  find  (weak)  solutions  ■u(tj)  to  (1)  for  1  <  i  <  n  <  oo. 
From  the  solutions  {u(tj,  z)}"=1  we  can  compute  the  average  temperature  change  between 
consecutive  time  steps  at  the  heat  sink  interface  (boundary  y2).  We  define  the  average 
temperature  change  from  tj_x  to  t*  by 

T-t  =  f  Tr2  ( u(U )  -  u(ti-i))(z)  dS2  (28) 

I  /  2  I  -7  72 

for  7  =  2,...,  n.  Recall  Tr2  is  the  continuous  trace  operator  from  iT1  (0)  — >■  H %  (y2)  defined 
by  (Tr2  /)(^)  =  /(z)|72.  We  choose  to  look  at  the  data  this  way  in  order  to  relate  our  model 
data  to  our  experimental  data. 

Before  using  the  experimental  data  we  have  collected  in  our  parameter  estimation  prob¬ 
lem,  we  first  used  generated  data  for  proof  of  concept.  We  generated  data  by  solving  (1)  at 
the  eight  times  corresponding  to  our  data  with  the  parameter  values  from  Table  1.  We  call 
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these  solutions  u(ti,z),  n(l2.  : );  . . . ,  u(t8,z).  We  define  the  vector  of  average  temperature 
changes  T  =  [T2  T3  •  •  •  X8]  where 

Tj  =  |^j  J  Tr2  (u(tj)  -  uitj-iMz)  dS2 

for  2  <  j  <  8.  Furthermore,  we  define  T(q)  to  be  the  vector  generated  by  the  solution  using 
the  specified  parameter  (or  parameters)  q.  which  we  will  denote  u(t,z;q),  and  (28),  i.e. , 

Tj{q )  =  [  Tr2  (■ u{tj ;  q)  -  »(/,  ?))(*)  d<S2 

1 72 1  J72 

for  2  <  j  <  8  and  T(r/)  =  [T2(g)  T3(q)  ■  ■  ■  T8(q)\  .  We  will  assume  any  unspecified  parameters 
are  given  by  the  values  in  Table  1. 

We  first  tried  to  estimate  a  constant  value  for  the  thermal  conductivity  of  the  aluminum 
particles  (for  examples  with  uniformly  distributed  particles  of  equal  size)  that  best  matched 
our  generated  data.  This  parameter  estimation  problem  was  unsuccessful  even  in  this  sim¬ 
plest  of  cases.  In  comparison,  we  were  successful  in  estimating  the  constant  value  for  the 
thermal  conductivity  of  the  silicone  polymer  that  best  matched  our  generated  data.  These 
results  were  not  surprising,  however,  when  the  graphs  of  the  cost  functions  for  the  parameter 
identification  were  plotted  as  a  function  of  the  parameter.  The  cost  function  for  the  particle 
parameter  was  jagged  with  no  clear  minimum,  whereas  the  cost  function  for  the  silicone 
parameter  was  smooth  with  a  clearly  defined  minimum.  For  further  information  on  our 
constant  parameter  estimation  problem  see  [8]  and  for  general  inverse  problems  see  [3,  4,  5]. 

4.2  Estimating  the  Thermal  Conductivity  Parameters  as  Random 
Variables 

The  theoretical  basis  of  this  approach  can  be  found  in  [1]  (see  also  [9,  12]).  While  this 
formulation  allows  us  to  estimate  the  distributions  for  both  of  the  thermal  conductivity 
parameters,  kp  and  ks ,  at  the  same  time,  we  will  estimate  them  individually  here.  We  first 
assume  the  parameter  kp  is  a  realization  for  a  normally  distributed  random  variable  with 
mean  jtp  and  variance  ap  and  attempt  to  estimate  its  distribution.  For  now,  we  will  hold 
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all  other  parameters  constant.  In  order  to  ensure  all  values  in  the  distribution  for  kp  are 
positive  we  will  use  a  truncated  Gaussian  distribution  and  require  nP  —  3 ap  >  0.  Thus  the 
probability  density  function  for  the  random  variable  is  given  by 


f{x)  = - =  (29) 

0.9974  op\/7h{ 

for  x  e  [np  —  3 ap,  fi p  +  3crp\.  If  we  define  qp  =  (i fip ,  crp),  we  can  define  the  expected  value  of 
T  by 


E[T[U,qp)\Pp  = 


r*  fJjp-\-3(Jp 


flp  S(Tp 


T(U ;  x)  1  ^  d 

v  ’  0.9974  <7„V2vr 


where  Pp  is  the  probability  distribution  function  arising  from  (29)  with  mean  pp  and  variance 
a?.  Thus  the  “best  fit”  parameter  q*  is  the  solution  to  the  least  squares  problem 


min  J{Pp,  T)  =  min  J(Pp )  =  min  V  |£[T(i/:;  qp)\Pp\  -  T,\ 
qP£Q  qP£Q  qP£Q 

i=i 


where  Q  =  M+  xl+  with  the  additional  restriction  pp  —  3 ap  >  0. 

For  proof  of  concept,  we  used  the  same  generated  data  described  in  Section  4.1.  While 
we  realize  the  data  was  not  generated  with  kp  normally  distributed,  we  would  consider  a 
successful  parameter  estimation  to  have  the  mean  of  the  distribution  near  the  true  value  for 
kp  and  small  standard  deviation.  We  carried  out  numerous  estimation  trials  and  in  Table  2 
we  present  values  for  two  of  these  minimizations  of  (30).  We  used  Matlab’s  constrained  mini¬ 
mization  routine  fmincon.  The  function  fmincon  uses  a  Sequential  Quadratic  Programming 
method.  The  three  main  steps  of  this  algorithm  are  the  solving  of  a  Quadratic  Program¬ 
ming  subproblem,  the  updating  of  the  Hessian  matrix  of  the  Lagrangian  solution,  and  the 
calculation  of  a  merit  function  and  line  search.  A  complete  description  of  this  method  can 
be  found  in  [15]. 

The  initial  guess  is  denoted  by  qPo  =  {nPo,  <jPo  )  and  the  best  parameter  fit  found  by  Matlab 
is  denoted  by  q*  =  {p*,a*).  For  all  of  these  minimizations,  n  =  5  in  (30). 

In  all  our  tests,  the  mean  of  the  “best  fit”  distribution  did  not  closely  resemble  the  actual 
parameter  value,  kp  =  217  and  the  standard  deviation  was  not  small.  In  all  cases,  the  optimal 
parameters  are  not  far  from  the  initial  guesses.  While  we  used  a  relatively  small  number 
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fho 

aPo 

J{(lpo ) 

k*P 

ap 

WP) 

1000 

100 

0.0197193 

1242.66 

297.228 

0.0153633 

250 

50 

0.0391095 

250.001 

50.0005 

0.0391080 

Table  2:  Estimating  distribution  of  kp 


for  n  and  a  constant  distribution,  we  would  still  expect  better  results  if  the  inverse  problem 
were  well  behaved. 

In  contrast,  suppose  we  instead  view  ks  as  a  realization  for  a  normally  distributed  random 
variable  and  estimate  and  as  for  the  distribution.  As  before  we  used  a  truncated  Gaussian 
to  ensure  all  possible  values  are  positive,  with  probability  density  function 


f(x)  = 


d-(x-/j,s)2 /2a2 


0.9974  oss/^ 

for  x  G  [//,,  —  3 as,  jis  +  3 cra\,  fis  —  3 as  >  0.  If  we  define  qs  =  ( /is. ,  as ),  we  can  define 

p  /^s+3<r,s 

mu-,  =  /  m;  xIsstiGb  dx 

j  Hs-ZcTs 


(31) 


where  Ps  is  the  probability  distribution  function  arising  from  (31)  with  mean  //,,  and  variance 
a2.  Thus  the  “best  fit”  parameter  q*  is  the  solution  to  the  least  squares  problem 


n 

min  J(PS,T)  =  min  J(PS)  =  min  V]  |£[T(tg  qs)\Ps]  -T)|2  (32) 

qs£Q  Qs&Q  qs&Q  A ' 

i=l 

where  Q  =  M+  xl+  with  the  additional  restriction  //p  —  3 ap  >  0. 

In  Table  3  we  see  values  for  one  minimization  of  (32)  (again  we  carried  out  multiple 
tests,  obtaining  similar  results).  We  used  Matlab  and  the  generated  data  as  before.  The 
initial  guess  is  denoted  by  qSo  =  (juSo,aSo)  and  the  best  parameter  fit  found  by  Matlab  by 
q*  =  (/U*,a*).  For  all  of  these  minimizations,  n  =  5  in  (32). 

In  this  example,  note  the  mean  [i*  of  the  distribution  is  very  close  to  the  actual  value 
of  ks  =  0.12  W/mK  used  to  generate  the  data  and  the  standard  deviation  is  small.  This 
indicates  that  despite  the  fact  we  are  using  a  small  number  for  n,  the  parameter  estimation 
algorithm  performs  well  when  estimating  ks  (or  its  distribution). 
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Tso 

as  0 

J(Qs0) 

T*s 

as 

m) 

1 

0.25 

19.294 

0.12036 

2.5087e-05 

1.4140e-02 

Table  3:  Estimating  distribution  of  ks 

5  Sensitivity  Equations 

In  the  previous  section  we  found  that  small  changes  in  the  parameter  for  the  thermal  con¬ 
ductivity  of  the  particles,  kp ,  seem  to  have  little  impact  on  the  solution  u(t,z).  In  contrast, 
small  changes  in  the  thermal  conductivity  of  the  silicone,  ks,  appear  to  significantly  change 
the  solution  u(t,z).  Sensitivity  equations  allow  us  to  anticipate  and  quantify  how  changes 
in  the  parameters  affect  changes  in  the  solutions.  Thus,  we  turn  to  sensitivity  equations  to 
investigate  whether  kp  and  ks  are,  in  fact,  influencing  the  solutions  as  we  suspect.  Readers 
are  referred  [19,  20,  21,  22]  for  more  information  on  sensitivity  equation  methods  and  their 
use  in  inverse  problem  methodology. 

5.1  An  Abstract  Derivation  in  Terms  of  Frechet  Derivatives 

In  order  to  derive  the  sensitivity  equations,  we  must  formally  differentiate  our  system  of  equa¬ 
tions,  including  the  boundary  conditions,  and  then  interchange  the  order  of  differentiation. 
Before  explicitly  following  this  procedure,  we  first  present  a  framework  that  rigorously  justi¬ 
fies  the  derivation.  This  derivation  relies  on  the  Implicit  Function  Theorem  and  a  corollary 
to  the  Implicit  Function  Theorem,  both  of  which  we  state  here. 

Theorem  5.1  (Implicit  Function  Theorem)  [7,  Theorem  3.1.10,  p.  115]  Let  X ,  Y,  and 

Z  be  Banach  spaces.  Suppose  f(x,  y )  is  a  continuous  mapping  of  a  neighborhood  U  of  (x0,  y0 ) 
in  X  x  Y  into  Z,  f{xo,yo)  =  0  and  fy{xo,yo)  exists,  is  continuous  in  x,  and  is  a  linear 
homeomorphism  of  Y  onto  Z.  Then  there  is  a  unique  continuous  mapping  g  defined  in  a 
neighborhood  U\of  x0,  g  :  U\  — >■  Y ,  such  that  g(x0 )  =  y0  and  f(x,g(x))  =  0  for  x  G  U\. 

Corollary  5.2  [7,  Corollary  3.1.11,  p.115]  If,  in  addition  to  the  hypothesis  of  the  Implicit 
Function  Theorem,  fx(x,y )  exists  and  is  continuous  for  ( x,y )  near  (x0,y0),  then  the  function 
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g(x)  is  continuously  differentiable  for  x  G  U\  and 


g'(x)  =  -[/„(*, j(z))]  ‘/AjM). 


(33) 


A  proof  of  Theorem  5.1  and  Corollary  5.2  can  be  found  in  [7]. 

In  Section  3.2  we  found  there  exists  a  unique  weak  solution  to 


u  +  Au  =  F 
u(0)  =  u0 


(34) 


where  F  is  defined  in  (9)  and  A  is  given  by  (8).  Without  loss  of  generality,  we  assume  u ( 0 )  = 
0.  (If  not  use  a  simple  change  of  variables,  u  =  u  —  u(0).)  We  are  interested  in  examining 
the  sensitivity  of  these  solutions  with  respect  to  the  thermal  conductivity  parameter  k,  or, 
more  specifically,  with  respect  to  parameters  q  in  a  parameterization  k(q)  of  the  thermal 
conductivity.  Let  Q  be  the  space  of  all  possible  parameter  values,  y  =  L2(0,T;V*)  and 
IX  =  L2(0,T;  V)  fl  i71(0 ,  T;  V*).  The  norm  on  IX  is  given  by 


I  I  1 1 2  mm2  |i||2 

I  l$Jlu  =  IIC  li2(0,T;V)  +  I  M  l//1(0,T;V*) 

for  v  G  IX  (see  [13,  p.  102]). 

Note  F  as  defined  in  (9)  does  not  depend  on  q.  The  thermal  conductivity  does  depend 
on  q ,  k  =  k(q)  :  Q  — >■  R  and  so  by  (8)  we  have  A  =  A(q)  :  Q  — )■  2z?(V,  V*),  with 

{Am,  i’)v*y  =  <-*(9 W,  VV>):k  +  h  I  ( Tr2  d)(z)(Tr2  f)(z)  dS2 
9  J  72 

for  (f).  'W  G  V.  Thus  for  each  q  G  Q  there  is  an  associated  weak  solution  u(  ■ ;  q)  G  IX. 

Define  M  :  Q  — >■  Jf(U,  y)  by 


[M(q)]v  =  v  +  A(q)v. 

Since  (34)  possesses  a  unique  solution  for  each  F  G  I2(0,T;V*),  we  see  that  M (q)  is  invertible 
for  each  q  G  Q.  Note  also  M(q)  is  linear  in  f,  i.e. ,  [M(q)](vi  +  v2)  =  [M(q)\v i  +  [M(q)\v 2. 
Define  N  :  Q  — )■  y  by 

N(q)  =  F. 
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Finally,  define  G  :  'll  x  Q  — >■  y  by 

G(u,q)  =  {[M{q)]u)  -  [iV(g)], 

and  associate  u(-;q )  with  the  pair  (u({%q),q)  satisfying  G(u(-  \q),q )  =  0  in  the  y  = 
L2(0,T;V*)  sense.  Note  G  induces  a  natural  mapping  from  Q  to  'll  given  by  q  ^  u(  -;q). 
Also,  u  G(u(-,  q),  q)  is  an  affine  map  from  IX  to  ^ . 

Fix  q0  G  Q  and  let  u o  =  uo(-;qo)  G  'Ll.  In  order  to  characterize  the  sensitivity  at  q o 
we  need  the  operator  Dqu(q0 )  G  Jzf(Q,U).  We  assume  (u0,qo)  =  («0(  • ;  qo) ,Qo)  satisfies 
G(uo,  qo)  =  0  in  the  y  =  L2(0, T;  V*)  sense. 

Lemma  5.3  The  partial  Frechet  derivative  duG(uo,qo)  ■  U  — >■  y  exists  and  is  given  by 
duG{u0,q0 )  =  M{q0)  G  Jf(U,y). 

Proof:  For  h  /  0  G  U, 

\\G(uq  +  h,  q0)  -  G{u0 ,  q0)  ~  M(q0)h\\y 

=  || [M(qo)](uo  +  h)~  [N(q0)]  -  ([M(9o)]«o  -  [N(q0)])  ~  [M(g0)]/i|  |y 
=  \\[M{qo)]h-[M{qo)}h\U 
=  0. 

We  assume  the  function  k  :  Q  — >■  M  is  Frechet  differentiable  at  qo  and  so  we  define 
Q^J?(Q,J?(U,y))  by 

[^£(q0)u\  =  g/{q0)u 

where  .</  :  Q  — >■  Jzf  (Q,  Jz?(V,  V*))  is  given  by 

(.W((/o)<A0)v*,v  =  {^Dqk{qo)V(j),V4’)si 

for  0,0  G  V.  Note  that  X  :  Q  — >■  M  and  so  Dqk(q0)  G  Jzf(Q,R). 

Let  W  =  2z?(U,  ^).  Since  ^  and  'll  are  Banach  spaces,  is  a  Banach  space  as  well. 
Thus  we  can  define  a  norm  on  W  by 

imi4t-'=  sup  1 1 Tv  1 1  y 

Mlu=i 
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for  T  G  .  Recall  for  (eV: 


V*  =  sup  ((,  '0}v*,V 

1 1^1 1  V=i 


and  for  <b  G  V 


M\%  =  l|V0||?c  +  Will,. 


Recall  from  (2)  that  g  is  positive  and  finite  and  note  that  for  any  <p  G  V, 


IIWII  12(q)  ~  (V</>,  V0)L2(n) 

=  -gV(f>,  V0)L2(n) 

=  (-V0,  Vd>)^ 

9 

<r~l  \m\i< 

<RlM\% 

and  so  1 1 V<f>\ \L2(n]  <  sJr^  ||0||y 

Lemma  5.4  M  is  Frechet  differentiable  at  qo,  with  [.Jf(qo)ti]h  =  [DqM(q0)h]u.  Further¬ 
more,  ^(qf)  is  a  bounded  linear  operator. 


Proof:  For  /i/OeQ  we  want  to  show 

lim  — l—  \\M{q0  +  h)  -  M(q0 )  -  [^(qo)]h\\&  =  0. 

h-^°  Iloilo 

By  definition  of  the  norm  in  FF . 

\\M{q0  +  h)  -  M{q0)  -  [^{q0)]h\\&  =  sup  ||[M(g0  +  h)u  -  M(q0)u  -  [,Ff{q0)u]h\ ||y. 

IMIu=i 

For  any  u  G  R  with  ||u||u  =  1, 

1 1  [M  ( g0  +  h)  u  -  M  ( q0)u  -  [Jt  ( q0 )  u]  h\  \  \  \ 
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Thus 


Jo 


\\{[M{q0  +  h)]u){t)  -  ([M{q0)]u){t)  -  {[^(q0)u\h)(t) ||^  dt 


| u(t)  +  A(q0  +  h)u(t)  -  [u(t)  +  A(q0)u(t)]  -  srf  (q0)u(t)h\\%~  dt 


J  o 


sup  (A{q0  +  h)u{t)  -  A{q0)u{t)  -  srf ( q0)u{t)h ,  V’)v-,v  dt 

'0  Il\’=1 

T  1  1  1 

sup  (-k(q0  +  h)\7u(t ) - k(qo)\7u(t) - Dqk(q0)h'Vu{t ),  dt 


Jo 


9 


9 


9 


< 


Jo 


sup  ||fc(5o  +  /i)  -  &(9o)  -  DgkiqojhWlc  { Vu{t ),  V0)| .»(n)  dt 


l7  — 1 


<  sup  ||fc(g0  +  ft)  -  k(q0)  -  ||V0|||2(n)  /  ||Vu(t)|||2(n)  dt 

||'0||v  =  l  JO 


<  sup  ||fc(5o  +  h)  —  k(qo)  —  Dgk(qo)h\\oc  RL  ||0||v  /  i?L  ||u(t)||v  dt 
1101k—1 


—  \\k{qo  +  h)  k{q0)  Dqk(qo)h\\00  RL  \ \u\ |.£2.(0iT;v) 
<  \\k{q0  +  h)  -  k{q0)  -  Duk{q())h[  'i  R IMIu 


\\M(q0  +  h)~  M{q0)  -  ■ //(qA'Ii  '  <  \\k{q0  +  h)  -  k{q0)  -  Dqk{q0)h\\oo  RLl . 


We  know 


lim  — - 

/i~>o  h  q 


||&(go  +  h)  ~  k((lo)  ~  Dgk^hWoo 


since  we  assume  k  is  Frechet  differentiable.  Thus 


0 


lim  — ||M(f/0  +  h)  -  M(q0)  -  V/(qo)\h  j:  =  0 

h^0\\h\\Q 

and  so  M  is  Frechet  differentiable  with  DqM  =  ,M  and  €  Jzf(Q,  Jf(U,  y )) 

Af(Q,  X). 
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Lemma  5.5  If  DqN(q0)  :  Q  — >■  y  and  DqM(q0 )  :  Q  — >■  Ff(U,y)  exist  in  the  Frechet  sense , 
then  the  partial  derivative  of  G  with  respect  to  Q  at  (uo.qo)  exists  and 
[dqG{u0,q0)]  £  Ft f(Q,y )  is  given  by 

[dqG(ua,  q0)]  =  Jt{ q0)u0  -  DqN{q0). 

Proof:  By  Lemma  5.4.  DqM(an)  =  ^{qo)  £  Jzf(Q,  J?(U,y)).  Since  N(q0)  has  no  dependence 
on  q,  DqN(q0 )  is  the  zero  operator  and  so  DqN(q0 )  £  Jf(Q,y).  Thus  ^(q0)u0  —  DqN(q0)  £ 

Ff(Q,y). 

For  h  ^  0  £  Q, 

\\G{u0,qo  +  h)—G(u0,q0)  —  (q0)u0  —  DqN(q0)]h\\y 

=  || [M(q0  +  h)]u0  -  N(q0  +  h)  -  ([M(g0)]«o  -  N(q0)]) 

-  ([ DqM{q0)h]u0  -  DqN(q0)h\\v 

<  || M{q0  +  h)  -  M(q0)  -  DqM(q0)h\\^\\u0\\u 
+  \\N(q0  +  h)-N(q0)-DqN(q0)h\\y 

Since  DqM(q0 )  exists  in  the  Frechet  sense, 

lim  jjTjj— | |M(g0  +  h)  -  M{q0 )  -  DqM(q0)h\\$-  =  0. 
h^o\\h\\Q 

Moreover,  since  N  has  no  dependence  on  qi 


N(q0  +  h)  -  N(q0 )  -  DqN(q0)h  =  0. 

Thus 

lim  — —  \\G{uo,qo  +  h)  -  G{u0,q0)  -  q0)u0  -  DqN{q0)]h\ |y  =  0 

h^°\\h\\Q 

and  so  [dqG{u0,  q0)]  =  Jt{ q0)u0  ~  DqN{q0). 


Theorem  5.6  Let  Q0  be  a  subset  of  the  interior  of  Q  and  let  1X0  be  a  subset  of  the  inte¬ 
rior  of  U.  Fix  qo  £  Qo  C  Q.  Suppose  there  exists  a  unique  uq(  ■ ;  qo)  £  Uo  C  U  such  that 
G(u0(  •;<7o)i(7o)  =  0-  //[M(g0)]_1  exists  in  Jzf(y,U)  and  if  DqM(q0)  =  -#(<7o)  and  DqN(q0 ) 
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exist  in  the  Frechet  sense  in  SY{Q,  JT)  and  Jf(Q,y)  respectively,  then  the  sensitivity  opera¬ 
tor,  s  =  Dqu(  ■ ;  qo)  E  F?(Q,U),  exists  and  satisfies 

M(q0)s  =  ~^{q0)u0  +  DqN(q0) 

Proof:  In  the  Implicit  Function  Theorem  5.1,  we  set  X  =  Q,  Y  =  It,  and  Z  =  y.  The  function 
f{x,y )  in  Theorem  5.1  is  our  G(u,q ),  and  by  assumption  G(uo,qo )  =  0.  By  Lemma  5.3,  we 
know  duG(u0,  q0)  £  Jz f(U,y).  Since  M(q)  is  invertible  for  each  q  £  Q.  [M(q0)]~L  £  Jf(y,U). 
Thus  u(  ■ ,  q)  satisfies  G(u(  ■ ,  q),  q)  =  0  for  q  £  Q0. 

By  Lemma  5.5,  we  know  dqG(u0,q0 )  £  Jf(Q,y).  Furthermore,  by  Lemma  5.4  we  know 
DqM(q0 )  =  and  so  by  Corollary  5.2,  M(q0)s  =  —^(qo)uo  +  DqN(q0).  In  fact,  since 

DqN(q0)  is  the  zero  operator,  M(q0)s  =  —^(q0)uo- 


5.2  Sensitivity  to  the  Particle  Thermal  Conductivity 

5.2.1  Derivation  of  Sensitivity  Equations 

In  Section  3.1  we  established  the  system  (1)  has  a  weak  solution  u(t,z).  If  we  assume  the 
thermal  conductivity  is  dependent  on  the  parameter  q ,  then  q  i— »  k(z;  q)  is  given  by 


k(z \q) 


(35) 


We  want  to  consider  the  thermal  properties  of  the  composite  silicone  as  we  hold  the  constant 
qs  fixed  and  let  the  constant  qp  vary  over  a  range  of  admissible  material  values.  Thus  any 
weak  solution  of  (1)  will  have  the  form  u(t,  z;  q).  By  studying  the  sensitivity  of  the  solution 
u  to  changes  in  qp.  we  can  study  the  sensitivity  of  the  system  to  the  thermal  conductivity  of 
the  particles. 

In  order  to  derive  the  sensitivity  equations,  we  will  formally  differentiate  (1)  with  respect 
to  qr>.  interchange  the  order  of  differentiation  and  define  the  sensitivity.  Our  analysis  in 
Section  5.1  guarantees  the  existence  of  these  derivatives  and  that  the  resulting  sensitivity 
equation  can  be  rigorously  interpreted  in  terms  of  an  associated  weak  or  variational  system. 
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First,  formally  differentiating  g(z)ii(t,z )  =  V  •  (k(z)Vu(t,  z))  with  respect  to  qp  yields 
9{z)-^~{u{t,z-q))  =  A(v  •  (k(z;q)Vu(t,#fq))). 

UQp  VQp 

Switching  the  order  of  differentiation  we  see 

c)  r)k  i 

g(z)-^;{-^-(t,z;q))  =  v  •  {—{z;q)Wu{t,z;q))  +  V  •  (k(z;  q)V—(t,  z;  q)). 
at  oqp  oqp  oqp 

If  we  then  define  the  sensitivity  to  qp  as  s(t,  z\  q)  =  J^(£,  z;  q)  and  substitute  into  the  previous 
equation  we  have 

d  k 

g{z)s(t,  z;  q)  =  V  •  (  — ( z;  q)Wu{t,  z;  q ))  +  V  •  {k{z]  q)Vs{t,  z;  q)). 

oqp 

It  is  also  necessary  to  differentiate  the  boundary  conditions  with  respect  to  qp.  To  do 
this,  note 


d_ 

dqp 


{k{z;  q) 


du 

dn 


(t,z\q)) 


dk  .  .du  d  du 

w— (a  < i )  — —  ( / .  q)  +  k(z:  '/ItHtt— (/•  A  q)> 


dq, 


dn 


d  n  dqp 


,  .  .  ds  dk  sdu 

k{z]  q)  —  {t,  $]  q)  +  w—  (2;  q)-^-(t,  z ;  q). 


dn 


dqp 


dn 


We  will  assume  the  source  flux  So  and  the  initial  condition  T  are  independent  of  qp. 
It  is  important  to  note  that 


dk 

dqp 


( z\q ) 


0  2  g  Q, 

1  z  G  fip 


and  so  j^-(z;q)  G  L°°(Q).  Furthermore,  J^-( z;q )  =  0  since  we  are  holding  qs  =  qs  fixed. 
Since  q  =  (, qp,qs ),  Dgk(q)  =  (|^,|J)  =  (|^,0).  Also,  since  u  G  L2(0,T;  V),  we  know 
Vu  G  L2(0,  T;  L2(Q)2).  Then  if  we  define 

d  k 

fit,  A  (l)  =  w—  {z;  q)Vu{t,  z ;  q),  (36) 

dqp 

we  have  /(•,•;  9)  G  L2(0,  T;  L2(kl)2). 
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Thus  we  formally  have  the  following  system  for  our  sensitivity  equation: 


g{z)s{t,,  z;  g)  =  V  •  ( k(z ;  g)Vs(£,  z;  g))  +  V  •  f{t ,  z;  g) 

^w®mw)|74  =  -J^;9)£(M;9)l74 
<  £(«w)fKM;g)|72  =  +  MMw))|7*  ^ 

^w)f^(Mw)|7l  =  ^U:'/)^(/.  ?W)  7l 

k(z:(j)^(l.r.:(j)"/):i  =  ~^{z\q)^{t,  Z]  q)  |7s 
s(0,  2;  g)  =  0 

Note  the  solution  u(t,  2:  g)  to  (1)  acts  as  part  of  the  forcing  term  /  on  the  solution  to  (37). 

5.2.2  Weak  Solution  to  Sensitivity  Equation 

We  refer  to  the  spaces  H,  V,  IK,  and  V  defined  in  Section  3.1.1.  Define  a  sesquilinear  form 
P  :  V  x  V  -»  R  by 

VO  =  <-*W,  +  h  I  (Tr2  P)(z)(Tr2  4>){z)  dS2  (38) 

9  J  72 

and  note  /3  defines  an  operator  01  G  Jzf(V,  V*)  where  (010,  V’)v*,v  =  P(<P,  V’)-  We  observe  that 
/3  =  /3(?)  is  just  the  parameter  dependent  sesquilinear  form  a  of  (8)  and  01  =  01(g)  is  the 
analog  of  A  generated  by  a. 

Define  T :  [0,T]  ->■  V*  by 

=  (39) 

where  /  is  given  in  (36).  We  shall  see  shortly  that  01  and  £F  are  the  operators  we  need  to 
find  the  weak  or  variational  form  of  (37). 

Suppose  s  solves 

s  +  01s  =  T  in  V* , 

i.e. ,  for  all  ip  e  V, 

(s(f)  +  0ls(t)  -  T(£),  ip)v*y  =  0. 
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By  definition 


(s{t)  +As(t),4>)v.y  =  (T(t),V’)v*,v 


and  thus 


(s{t),  4')v*y  ={ — kVs{t),  V'0)jc  -  h  /  {Tr2  s(t))(z)(Tr2  t/’)(z)  dS2 
9  .7  72 

+  <-7/(«),  w)M. 


(40) 


Note  that  the  result  in  Theorem  5.6  says  the  sensitivity  operator,  s,  satisfies  M(q)s  = 
— (<?)«,  i.e. , 

s(t)  +  A(q)s(t)  =  -g/(q)u, 

which  by  definition  is 

(s{t),  V’)v*,v+(-&Vs(i),  VV>)rK  +  h  f  {Tr2  s(t))(z)(Tr2  4>){z)  dS2 
9  j  72 

Clearly  this  is  equivalent  to  (40)  and  hence  our  formal  and  rigorous  derivations  result  in  the 
same  system. 

Returning  to  equation  (40),  we  see  it  is  equivalent  to 

<s(t),V’)v*,v  ={-kVs(t),  Vt/’)z’  -  h  [  (: Tr2  s(t))(z)(Tr2  4>)(z)  dS2 


’  72 


(  ./(/).  V  r)/;.! 


by  the  definition  of  the  fit-norm. 

Now,  if  s  G  L2(0,T;V),  /,;  Vs  G  L2(0,T;V),  and  /  G  L2(0,T;R),  using  the  Divergence 
Theorem  and  the  identity 

V  •  {k{z)Vs{t,,z)ilj(z))  =  (V  •  (k(z)Vs(t,  z))i/’(z)  +  k(z)(Vs(t,  z)  •  Vil’(z)) 


in  the  previous  equation,  we  see 

<-(/)• R)v.v  =(V  •  (fcVs(t)),y'-)i2  -  /  (Tr0  (n  •  fcVs(*)))(*r)(7Vn  </’)(«)  ^ 

Jan 
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-  h  (Tr2  s(t))(z)(Tr2  4>){z)  clS2 

J  72 

-  [  (' Trn  (n  ■  f(t)))(z)(Trn  p>){z)  clS  +  (V  •  /(£), 

Jan 

where  T/’n  is  the  continuous  trace  operator  mapping  /  E  H1(Q )  —>■  H^(dQ)  with  (Tr q  f)(z)  = 

f(z)\m- 

If  we  have  enough  smoothness  on  s  with  respect  to  the  time  derivative  (i.e. ,  s  E  T2 (0,  T;  IK) ) 


{s(t),4>)vy  =  (s(t),V’)rK  =  {gs{t),4')L  2 


and  so 


(gs(t)  V  •  (/.'V //(/))  V  •  /(/).  = 

/•  07. 

-  /  (Trn  n  •  (feVs(i)  -  —  Vu(t)))(T)(Trn  </’)(*)  dS 
-  h  f  ( Tr2  s(t))(z)(Tr2  xp){z)  clS2 

J  72 

for  all  xp  E  V.  We  know  (41)  holds  for  all  xp  E  V  =  TT1  (0) ,  thus  for  all  xp  E  Hq(Q)  , 

(gs(t)  ~  V  ■  (k\7u(t))  -  X7  ■  f(t),4’)L2  =  0. 

Since  Hq(Q)  is  a  dense  subset  of  T2(fi),  we  know 

gs(t)  —  V  •  (kVu(t))  —  V  •  f(t)  =0 

in  the  L 2  sense.  Thus  we  can  rewrite  (41)  as 

r  or 

-  /  (Tm  n  ■  (kVs(t)  -  —Vu(t)))(z)(Trn  v)[z)  dS 
Jan  adp 

-  h  [  ( Tr2  s(t))(z)(Tr2  rii’){z)  dS2  =  0. 


We  know  dVt  =  71U72U73U74  .  Thus  for  xp  G  =  {<P  E  H x(fi)  :  </>|72,73,74  =  0}  C  i71(fi), 
by  (42)  we  see 

r  or 

-  /  (Tn  n  •  (fcVs(i)  +  —  Vu))(^)(Tn  t/’)(*)  dSi  =  0 

7  7!  VQp 
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and  so  k(z)^(t,  z) |7l  =  --^-(z)^(t,z) |7l.  A  similar  argument  shows 
k{z)^{t,z) |73  =  and  k{z)^{t,z) |74  =  -  JJWff  (Mk 

For  ft )  E  Hyi  =  {<p  E  H 1(fi)  :  0|7l,73,74  =  0}  C  using  (42)  we  see 


r  f)U 

-  /  (Tr2  n  •  (fcVs(t)  +  —  Vu))(^)(Tr2  </’)(*)  dS2 
I  dqp 


'72 


-  h  /  (Tr2  s(t))(z)(Tr2  ft’){z)  dS2 


'72 


0 


and  thus  z)^(t,  z)\l2  =  —(-jjjj-(t,z)^(t,x)  +  /is(t,  ^))|72.  Thus  if  s  satisfies  (40),  and  / 
and  s  have  sufficient  additional  regularity,  s  provides  a  strong  solution  to  (37). 

Thus  (40)  is  the  weak  or  variational  form  of  (37),  and  hence  any  solution  s  of  (40)  (if  it 
exists)  is  a  weak  solution  of  (37). 


5.2.3  Well-Posedness 


We  next  establish  existence  of  solutions  to  parabolic  systems  of  the  form 


s  +  As  =  T  in  V* 

S(0)  =  -Sq. 


(43) 


We  use  the  Gelfand  triple  V  <— >■  =  “K*  >  V*  as  in  Section  3.2.  Since  ft  defined  in  (38) 

is  equivalent  to  a  defined  in  (8),  we  know  ft  is  V-bounded  and  V-coercive  uniformly  in  q  E  Q. 
We  also  observe  that  the  forcing  term  T  defined  in  (39)  satisfies 


Tg  L2(0,T;W). 


(44) 


In  order  to  see  that  (44)  holds,  recall  /  G  L2(0,T;  L2(ftl)2).  Thus  T(t)  :  V  — >■  R  and  is  linear, 
i.e.,  T(t)  G  V*  and  so  ?G  L2(0,T;  F). 

Given  the  above  hypothesis,  the  weak  or  variational  form  of  the  system  (43)  is 


I  +  ft (s(t),ft>)  =  (T(t),V>) 

[«(0)  =  «0 

for  (/>  G  V  and  ( • ,  • )  is  ( • ,  •  )v*,v-  Note  the  system  (43)  is  the  same  as  (45). 


(45) 
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At  this  point,  existence,  uniqueness,  and  continuous  dependence  on  the  forcing  function 
follow  using  an  argument  analogous  to  the  arguments  in  Section  3.2.  Thus  the  sensitivity 
system  is  well-posed. 


5.2.4  Matlab  Solutions 


Using  Matlab’s  PDE  Toolbox,  we  are  able  to  numerically  solve  the  sensitivity  equation 
system.  Since  we  are  interested  in  the  temperature  of  the  composite  silicone  at  the  heat 
sink  interface,  i.e.,  the  72  boundary,  we  are  interested  in  the  sensitivity  at  the  heat  sink 
interface  as  a  function  of  qp.  Recall  that  in  our  model  we  assume  the  particles  never  touch 
any  boundary  of  the  composite  silicone.  (A  reasonable  assumption  based  on  the  properties 
of  the  silicone  polymer.)  Thus  §^-{z;  q) |7j  =  0  for  j  =  1,  2,  3,  4  and  we  can  reduce  (37)  to 


g{z)s{t,  z;  q)  =  V  •  {k{z;  q)\7s{t ,  z;  q))  +  V  •  f(t,  Z]  q) 

k{z\q)^{t,z-q)  |74  =  0 

k{z\  z\  q)  |72  =  -hs{t,  z;  q)  |72 

k{z;,q)^(t.  z:q)  ,,  =  0 

k(z:q)§^{t,z-,q)  |73  =  0 

s(0,  z;  q)  =  0 


(46) 


To  solve  (46),  we  first  fix  a  value  for  qp  and  solve  (1).  Then  with  that  solution  we  can  solve 
(46)  with  the  same  fixed  value  of  qp.  We  will  define  the  average  sensitivity  at  the  boundary 
72  for  a  fixed  qp  at  time  t,  by 

«2 {U\qp)  =  —  f  Tr2  s(t4;qp)(z)  clS2  (47) 

72  Jy2 


where  Tr2  is  again  the  continuous  trace  operator  from  Hl(VL)  -z  (j2)  defined  by 
(' Jr 2  f)(z)  =  f(z) \l2.  Repeating  this  process  for  different  values  of  qp  G  Q,  we  can  then 
define  the  relative  average  sensitivity  to  the  thermal  conductivity  of  the  particles  by 


•52 r{ti]  qp) 


S'2 (U)  Qp) 

maxgp g q ( s2 ( ;  qp)) 
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We  fix  ks  =  qs  =  0.12  W/mK.  The  remainder  of  the  parameters  are  given  in  Table  1. 

In  Figure  4  we  depict  the  plots  of  s^rituQp)  as  a  function  of  qp.  We  let  qp  =  100,110, 

. . . ,  1100  and  solve  at  each  of  the  times  t2,  t3,  t4,  and  £5  (similar  graphs  for  £6,£7,  and  £g  are 
found  in  [8]).  Note  that  as  a  function  of  qr>.  there  is  little  variation  in  the  relative  average 
sensitivity  along  the  heat  sink  interface  at  each  time  step.  In  contrast,  in  Section  5.3.2  we 
present  plots  of  the  relative  average  sensitivity  with  respect  to  the  thermal  conductivity  of 
the  silicone,  and  find  that  there  is  a  substantial  variation  in  the  sensitivity  as  a  function  of 
qs.  Thus,  based  on  our  results  in  this  section  and  in  Section  5.3.2,  we  conclude  the  solution 
u  to  (1)  is  not  very  sensitive  to  the  thermal  conductivity  of  the  particles. 


Figure  4:  Relative  average  sensitivity  at  72  at  times  £ 2,  £3,  £4,  and  £5  as  a  function  of  qp 
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5.3  Sensitivity  to  the  Silicone  Polymer  Thermal  Conductivity 

5.3.1  Derivation  of  Sensitivity  Equations 

We  again  begin  with  (1),  and  now  we  assume  the  thermal  conductivity  k(z)  is  given  by 


k{z;q) 


qs  z  €  Ss 

qp  z  G  Slp 


where  qp  is  a  constant,  and  qs  is  varied  in  a  range  of  admissible  values.  Any  weak  solution 
of  (1)  will  again  have  the  form  u(t,  z;  q).  In  order  to  derive  the  sensitivity  equations  we  will 
differentiate  (1)  with  respect  to  qs. 

The  differentiation  follows  in  a  manner  similar  to  the  differentiation  in  Section  5.2.1.  We 
will  define  the  sensitivity  to  qs  as  w(t,z;q)  =  j^(t,z;q).  We  will  again  assume  the  source 
flux  Sq  and  the  initial  condition  T  are  independent  of  qs. 

It  is  important  to  note  that 


dk 

dqs 


{z;q) 


1  z  G  Sls 
0  z  eSlp 


and  hence  -§jj-{z;  q)  G  L°°(Q).  Furthermore,  §^-{z;  q)  =  0  since  we  are  holding  qp  =  qp  fixed. 
Thus,  since  q  =  {qp,qs),  Dqk{q)  =  =  (0,  -§^).  Recall  u  G  L2(0,T;V),  and  hence 

V'U  G  L2(0,  T;  L2(Sl)2).  Then  if  we  define 

d  k 

/,(/,  q)  =  T— (  q)Vu(t.  z;  q) 
we  have  /,(•,  -:q)  G  T2(0,  T;  L2(Q)2). 
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Thus  we  formally  have  the  following  system  for  our  sensitivity  equation: 


g(z)w(t,z;q)  =  V  •  (k(z;  q)Vw(t,  z;  q))  +  V  •  fs{t,z;q ) 
h(ziq)^(t.z:q)  ,,  =  -  J|(A  9)ff  (Mw)k 

k(z-,q)^{t,z-q)\^  =  ~{j^{z]  q)j£{t,  z;  q)  +  hw(t,  z;  q))\l2 
<  (48) 

k{z-q)^{t,Z]q)  |7i  =  —§^{z-,q)^{t,ziq)\n 
k{z;q)^{t,z;q)\l3  =  —§^{z;  q)§%(t,  z;  f/)|7s 
w( 0,  z;  q)  =  0. 

s. 

Using  an  argument  analogous  to  the  arguments  in  Section  5.2.2  and  Section  5.2.3  it  can 
be  shown  there  exists  a  unique  weak  solution  w  to  (48)  and  that  the  problem  is  well-posed. 


5.3.2  Matlab  Solutions 

As  in  Section  5.2.4  we  can  use  Matlab  to  solve  the  sensitivity  equation  (48)  for  different 
values  of  qs.  Since  we  assume  the  particles  never  touch  any  boundary  of  the  composite 
silicone,  J^-|7j.  =  1  for  j  =  1,  2,  3,  4  and  therefore  we  can  reduce  (48)  to 

g(z)w(t,  z;  q)  =  V-  ( k(z ;  q)Ww{t,  z;q))  +  V  ■  fs{t ,  z;  q) 
k(mq)^{t,z;q)\14  =  -S0(t) 

k{z;  q)^{t,  z;  r;)|72  =  -(Mr<x>  -  u{t ,  z;  q))  +  hw{t,  z;  g))|72 
<  (49) 

k(z;q)^(t,z;q)\7i  =  0 
k(r.:q)^(l.z:q)  ,3  =  0 
io(0,  z;  q )  =  0. 

s. 

We  fix  kp  =  qp  =  217  W/mK.  The  remaining  parameters  are  given  in  Table  1.  In  order  to 
implement  the  boundary  condition  for  y2  in  the  PDE  Toolbox,  we  will  use  the  average  value 
of  u{t,  z;  q)  on  y2  for  u(t,  z;  q)  in  k{z;  q)^{t,  $;  q)  |72  =  -(MToo  -  u{t ,  -jg;  q))  +  hw{t ,  z;  q))  |72. 
As  in  Section  5.2.4,  we  define  the  average  sensitivity  at  the  boundary  y>  for  a  fixed  qs  at 
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time  t,  by 


w2{U]qs)  =  —  Tr2  w(ti-,qs)(z)  dS2  (50) 

I  / 2  I  J 72 

where  w(ti,  qs)  is  the  solution  to  (49)  at  t,  for  a  given  qs  and  fixed  qp.  Repeating  this  process 
for  different  values  of  qs  G  Q.  we  can  then  define  the  relative  average  sensitivity  to  the 
thermal  conductivity  of  the  silicone  by 

w2 (U,Qs) 

rV'iiQ.s)  /  /,  \\- 

ma XqseQ(w2{Ul  qs)) 


Figure  5:  Relative  average  sensitivity  at  72  at  times  t2,  1 3,  1 7,  and  tg  as  a  function  of  qs 

In  Figure  5  we  present  several  plots  of  «;2?.(b;;  qa)  as  a  function  of  qs.  We  solved  (49)  with 
qs  =  0.02,  0.02, . . . ,  1.02  at  each  of  the  times  t2,  t3, . . . ,  t8.  In  each  these  graphs  it  clear  that  as 
a  function  of  qs  there  is  significant  variation  along  the  heat  sink  interface.  Thus  we  conclude 
the  solution  u  to  (1)  exhibits  significant  sensitivity  to  changes  in  the  thermal  conductivity 
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of  the  silicone.  Also,  we  note  the  significant  difference  between  the  graphs  presented  here 
and  those  in  Section  5.2.4.  Based  on  these  results,  increasing  the  thermal  conductivity  of 
the  particles  should  not  result  in  much  improvement  of  the  overall  thermal  conductivity  of 
the  composite.  (We  remark  that  this  agrees  with  initial  experimental  findings.)  Clearly 
the  composite  is  more  sensitive  to  the  thermal  conductivity  of  the  base  material.  Thus  a 
better  means  of  improving  the  overall  thermal  conductivity  of  the  composite  is  to  increase 
the  thermal  conductivity  of  the  base  polymer  used  in  the  composite. 

6  Conclusions 

We  have  presented  analysis  of  a  two  dimensional  model  based  on  the  composite  silicone  and 
the  data  collection  process.  We  summarized  numerical  findings  using  Matlab’s  PDE  Toolbox 
for  our  two  dimensional  model.  Matlab’s  PDE  Toolbox  allowed  us  to  accommodate  the 
oscillatory  coefficients  in  a  variety  of  particle  geometries.  (In  [8]  results  for  various  particle 
geometries,  including  random  and  uniform  geometries,  were  presented  in  some  detail.)  It 
was  found  that  the  geometry  of  the  composite  silicone  has  a  significant  impact  on  the  heat 
flux  at  the  interface  between  the  heat  sink  and  the  composite.  In  this  paper  we  have  given 
theoretical  results  for  general  geometries  and  some  numerical  findings  for  the  special  case  of 
a  uniform  geometry. 

We  have  given  existence  and  uniqueness  theorems  based  on  our  two  dimensional  model, 
and  have  shown  the  model  depends  continuously  on  parameters,  as  well  as  the  initial  data 
and  forcing  function.  We  have  presented  a  formulation  and  numerical  results  for  two  dif¬ 
ferent  parameter  estimation  problems:  estimating  parameters  as  constants  and  estimating 
parameters  as  realizations  of  random  variables.  Estimating  parameters  as  realizations  of 
random  variables  used  a  probability  based  approach,  and  we  have  provided  a  careful  theo¬ 
retical  framework  for  this  approach  in  a  separate  reference  [1].  All  of  these  results  readily 
hold  for  the  corresponding  three  dimensional  model. 

We  carried  out  several  numerical  experiments.  In  each  of  our  parameter  estimation 
formulations  (using  simulated  data  for  proof  of  concept)  in  two  dimensions,  with  a  uniform 
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geometry,  we  were  able  to  accurately  estimate  the  thermal  conductivity  of  the  silicone,  but 
not  the  thermal  conductivity  of  the  aluminum  filler  particles. 

After  deriving  the  sensitivity  equations,  we  studied  the  sensitivity  of  model  solutions  to 
both  the  thermal  conductivity  of  the  particles  and  the  thermal  conductivity  of  the  silicone. 
Our  numerical  results  clearly  indicated  the  solution  is  significantly  more  sensitive  to  the 
thermal  conductivity  of  the  silicone  than  to  the  thermal  conductivity  of  the  particles.  This 
supported  our  results  from  the  parameter  estimation  and  experimental  findings  to  date. 
Thus,  in  order  to  significantly  increase  the  thermal  conductivity  of  the  composite  silicone  (or 
any  composite  adhesive),  we  suggest  it  is  best  to  work  at  increasing  the  thermal  conductivity 
of  the  base  silicone  (or  base  adhesive).  However,  we  believe  there  is  still  a  great  deal  to  learn 
about  thermally  conductive  adhesives  using  the  methodology  in  this  paper  with  variable 
particle  geometries. 

While  we  have  not  presented  the  results  in  this  paper,  we  note  that  as  an  alternative 
to  using  Matlab’s  PDE  Toolbox  one  could  use  the  mathematical  theory  of  homogenization. 
Homogenization  [2,  6,  10,  17,  18,  24]  combines  the  oscillatory  coefficients  into  an  “aver¬ 
age”  or  “effective”  thermal  conductivity.  See  [8]  and  the  above  cited  references  for  further 
information  on  this  approach. 
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